% propagation constant

function [p coef] = func_poles(nx, CONSTS, span, plot_data)

%---------------------------poles------------------------------------------
    nz_low = span(1); 
    nz_high = span(2);
    num_intervals = nz_high*7;
    nz_vec = linspace(nz_low, nz_high, num_intervals)';

    func = @(nz_vec)abs(func_delta_from_nz(nz_vec, nx, CONSTS));
    
    if (plot_data)
        func_for_plot = abs(func_delta_from_nz(nz_vec, nx, CONSTS));
        figure; plot(nz_vec, func_for_plot, 'b.-'); title('|det|');
    end
    
    for i = 1:size(nz_vec,1)
        x0 = fminsearch(func, nz_vec(i));
        x(i) = x0;
    end
    
    j=1;
    x_pref(1)=x(1);
    for ii = 2:size(x, 2)
        if(abs(x(ii)-x(ii-1))>1e-3)
            j = j+1;
            x_pref(j)=x(ii);
        end
    end
    
%     ll = 0;
%     for l = 1:size(p_m_temp, 2)
%         if(p_m_temp(l)>=nz_low && p_m_temp(l)<=nz_high)
%             ll = ll+1;
%             p_m_pref(ll) = p_m_temp(l);
%         end
%     end
    
    p = sort(x_pref)';
    
%----------------------------coef------------------------------------------
    debug_output = false;
    for ii=1:size(p,1)
        p1 = p(ii); % was 1
        mat = matrix_of_boundary_conditions(p1, nx, CONSTS);
        det(mat);
        rank(mat);
        K = mat(1:3, 1:3);
        det(K);
        rank(K);
        C2 = (1/1i)*det(K);
        f = - C2 * mat(1:3, 4);
        u = K^(-1) * f;
        B1 = u(1);
        B2 = u(2);
        C1 = u(3);
        %C2
        p(ii)=p1;
        epsilon = 1e-10;
        %coef = [B1; B2; real(C1); real(C2)];
        %coef = clear_small_values([B2; B1; C1; C2], epsilon);
        coef1 = [B1; B2; C1; C2];
        coef(:,ii)=coef1;
    end
    
    eq = clear_small_values(mat * coef, epsilon);
    
    if(debug_output)
        fprintf('equations for p1  %e \n', eq);
    end

    if(debug_output)
        
        E_z_inside_a  = mat(2, 1)*B1 + mat(2, 2)*B2;
        E_z_outside_a = mat(2, 3)*C1 + mat(2, 4)*C2;

        H_z_inside_a  = mat(1, 1)*B1 + mat(1, 2)*B2;
        H_z_outside_a = mat(1, 3)*C1 + mat(1, 4)*C2;

        E_phi_inside_a  = mat(3, 1)*B1 + mat(3, 2)*B2;
        E_phi_outside_a = mat(3, 3)*C1 + mat(3, 4)*C2;

        H_phi_inside_a  = mat(4, 1)*B1 + mat(4, 2)*B2;
        H_phi_outside_a = mat(4, 3)*C1 + mat(4, 4)*C2;

        fprintf('r=a: E_z_inside   = %e + %ei      E_z_outside   = %e + %ei \n',  ...
                               real(E_z_inside_a),  imag(E_z_inside_a), ...
                               real(E_z_outside_a), imag(E_z_outside_a));
        fprintf('r=a: H_z_inside   = %e + %ei       H_z_outside   = %e + %ei \n',  ...
                               real(H_z_inside_a),  imag(H_z_inside_a), ...
                               real(H_z_outside_a), imag(H_z_outside_a));
        fprintf('r=a: E_phi_inside = %e + %ei       E_phi_outside = %e + %ei \n',  ...
                               real(E_phi_inside_a),  imag(E_phi_inside_a), ...
                               real(E_phi_outside_a), imag(E_phi_outside_a));
        fprintf('r=a: H_phi_inside = %e + %ei      H_phi_outside = %e + %ei \n',  ...
                               real(H_phi_inside_a),  imag(H_phi_inside_a), ...
                               real(H_phi_outside_a), imag(H_phi_outside_a));
    end
    
end
    
function det_abs = det_abs(nz, nx, CONSTS)

    det_abs   = zeros(size(nz, 1), 1);
    for i=1:size(nz, 1)
        p_i = nz(i);
        mat = matrix_of_boundary_conditions(p_i, nx, CONSTS);
        %det_vec(i) = det(mat);
        det_abs(i) = abs(det(mat));
    end

end


